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ABSTRACT 

We present a comparative analysis of several methods, known as local Lagrangian 
approximations, which are aimed to the description of the nonlinear evolution of large- 
scale structure. We have investigated various aspects of these approximations, such as 
the evolution of a homogeneous ellipsoid, collapse time as a function of initial conditions, 
and asymptotic behavior. As one of the common features of the local approximations, 
we found that the calculated collapse time decreases asymptotically with the inverse 
of the initial shear. Using these approximations, we have computed the cosmological 
mass function, finding reasonable agreement with iV-body simulations and the Press- 
Schechter formula. 

Subject headings: cosmology: theory — gravitation — dark matter — large-scale struc- 
ture of universe 



1. INTRODUCTION 

Large-scale structures are believed to have formed from the gravitational amplification of 
primordial perturbations. At its first stages, the process of gravitational clustering can be investi- 
gated using linear perturbation theory. However, as the universe evolves, nonlinear concentrations 
of mass arise. Many structures we see today correspond to fluctuations several orders of magni- 
tude higher than the mean density of the universe; for example, clusters of galaxies have typically 
Pciuster/Puniverse ~ 10 2 — 10 3 . For larger scales this ratio decreases, approaching unity in the largest 
structures. 
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As there is no analytical treatment for the nonlinear regime, iV-body simulations are often 
resorted to. The numerical simulations had an enormous development in the last decade (see 
Bertschinger 1998 and references therein), being able to reproduce many features of the large 
scale structure. However they do not always provide a clear insight of the physics of nonlinear 
gravitational collapse. Moreover, they are usually very time-consuming, making it difficult to scan 
a large part of the parameter space of the cosmological models. 

For this reason, semi-analytical methods have been devised to tackle such a complex problem. 
The first approximation developed to study the nonlinear regime was introduced by Zel'dovich 
(1970). There are now various approximation schemes to analyze different aspects of non-linear 
clustering, including extensions of the Zel'dovich approximation (for a review see Sahni & Coles 
1995). Among them, the so-called local Lagrangian approximations have been introduced rather 
recently. The basic feature of these local approximations is that the kinematical parameters in each 
fluid element evolve independently of those of other elements. Thus the time evolution of a self 
gravitating fluid is replaced by a set of ordinary differential equations. This comes at the expense 
of losing information about the positions of each fluid element. Only local quantities, such as the 
density contrast, shear and expansion rate can be determined. 

Due to their handy applicability compared to the numerical simulations, as seen in the case 
of the widely used Zel'dovich approximation, they deserve a closer investigation. For some of 
these methods, certain aspects of their performance and applicability have already been discussed. 
However, to the authors' knowledge, no systematic comparison among them has ever been done. 
We consider it worthwhile to analyze them in a unified way in order to exploit general properties of 
these approximations, clarifying their similarities and differences. It is also important to compare 
their performance in some practical applications. In this paper, we discuss the following four 
approximations, in addition to the Zel'dovich approximation: the Local Tidal Approximation (Hui 
& Bertschinger 1996), the Deformation Tensor Approximation (Audit & Alimi 1996), the Complete 
Zel'dovich Approximation (Betancort-Rijo & Lopez-Corredoira 2000) and the Modified Zel'dovich 
Approximation (Reisenegger & Miralda-Escude 1995). All of them intend to be applicable to 
the highly non-linear regime. To the best of our knowledge, these comprise all existent local 
approximations in the literature, that are exact for planar, spherical, and cylindrical symmetries 
(except for the Zel'dovich approximation). 

The paper is organized as follows: In section 2, we briefly review various local approximations 
in a unified way. In section 3, these methods are applied to several cases. First, we discuss the 
homogeneously collapsing ellipsoid. We then study their behavior under general initial conditions. 
Finally, we apply some of these approximations to the calculation of the cosmological mass function. 
We sum up our results and present conclusions in section 4. In the appendices we present useful 
formulae for the calculation of the mass function together with fitting formulae for the collapse 
time in the approximations considered here. 



-3- 



2. LOCAL APPROXIMATIONS 

Throughout this paper we will only consider the case of cold dark matter (CDM), which is 
assumed to be collisionless, at least on large scales. This is well justified since 80 to 90% of the 
matter that clusters is composed by CDM (Turner 2000, Durrer & Novosyadlyj 2000). Furthermore, 
as long as the trajectories do not intersect, we can treat the CDM as a pressureless fluid. 

We will be working in a matter-dominated flat universe (the Einstein-de Sitter universe, 
hereafter EdS). Recent observational evidences are consistent with a zero curvature universe (De 
Bernardis et al. 2000; Hanany et al 2000). Even if we had a non-flat universe we would only require 
that the curvature be negligible in the scales of interest. The assumption of matter dominance may 
seem unrealistic since the observations indicate that the universe is now dominated by a repulsive 
homogeneous cosmological term (Perlmuter et al. 1998; Riess et al. 1998; Zehavi & Dekel 1999). 
However, the energy density of this term decays more slowly than the matter density. In the case of 
a cosmological constant we would have pa = const, whereas for matter we have pu ot a~ 3 , where 
a is the scale factor of the universe. Since most structures form at a time when p\ <C pM, we can 
safely ignore the effect of the cosmological term on the collapse process. 

The peculiar motions in the universe are much smaller than the speed of light. For pertur- 
bations on scales smaller than the Hubble radius, we can use the Newtonian approximation to 
describe the gravitational clustering. The basic equations for nonrelativistic pressureless matter in 
a perturbed EdS universe are the Euler, the continuity and the Poisson equations (Bertschinger 
1996): 

1 dvj _ _d^_ 

a dr a 2 % dx % ' 

^ + a(l + 5)9 = 0, (2) 

UT 

where 5 = (p — p) jp is the density contrast V{ = (dxi/dr) /a is the peculiar velocity, 6 = dv l jdx 1 
is the expansion, cf> is the peculiar gravitational potential, and the time variable r is related to the 
cosmic time t (also known as proper time) by dr = dt/a 2 . The comoving coordinate x,- L is given in 
terms of the position by x% = n/a. The left hand side of equation (1) is simply (<i 2 Xj/<iT 2 ) /a 2 , 
so that it looks like the usual Euler equation (apart form the factor a~ 2 ). In an Einstein-de Sitter 
background the scale factor is proportional to r -2 . We set a = t~ 2 such that AiiGa 2 p = 6r 2 = 6/a. 
The present value of the scale factor ao is fixed to be unity. 

The Lagrangian coordinates qi are often used instead of the position Xi in nonlinear analyses. 
In terms of the convective derivative d/dr = d/dr\ x + Vid/dxi is simply given by the time 
derivative at fixed q: d/dr = d/dr\ q . The Lagrangian coordinates are chosen to be the initial 
comoving positions: q^ = lim a ^o r i/o- 
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The Jacobian matrix of the transformation X\ — > % 

J « = §> < 4 > 

is known as the deformation tensor. The velocity gradient dvi/dx^ can be expressed in terms of 

J i j 9jS 

3a* " a dr • W 

The density is given by p(x,t) = pj, where J is the determinant of Jij. It is easy to see that the 
continuity equation (2) is solved exactly with 5 = J~ l — I. 

Differentiating equation (1) with respect to Xj we find 

J-i^l = -_ ^_ (6) 
i k dr 2 a 2 W W 

whose trace furnishes Raychaudhuri's equation 

J-/^ = -4nGa^(J- 1 -l)- (?) 

This is a local equation for Jij in the sense that it has no spatial derivatives, although it is not 
sufficient for determining the nine components of the deformation tensor. Usually this equation is 
written in terms of the kinematical parameters, 6, Oij (shear) and u>i (vorticity), defined by: 



-q^~ — g & $ij + a ij + ^ij a ij — a ji i ^ij ~ e ijk ^ — (8) 

If the initial conditions have no vorticity, then we have uj 1 = during all the evolution, as long as 
the trajectories do not intersect. Here we will consider only the case of vanishing vorticity. 

Equations (1) to (3) form a set of nonlinear partial differential equations. However, for certain 
specific configurations the time evolution of the deformation tensor behaves as if each space point 
evolves independently from the others. One might then expect that for more general situations the 
locality may hold, at least approximately, for these variables. Accordingly, several methods have 
been introduced which are known as local approximations. In their framework, the influence of the 
neighbors may enter only through the initial conditions. 

In addition to the solution of the continuity and Euler equations, the local approximations 
discussed here will replace the essentially nonlocal exact equation (3) either by some Ansatz inspired 
on equation (7), or by local evolution equations for the second derivative of the peculiar gravitational 
potential 4> (see also Kofman & Pogosian 1995 for a discussion). 

One of the basic features of local approximations is that the eigenvectors of the deformation 
tensor do not change with time. Thus, once diagonalized, remains diagonal in the same frame, 
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along all the evolution. This condition is either assumed from the beginning or appears as a 
consequence of the approximation introduced in the evolution equations. Actually, this assumption 
is not strictly consistent with the evolution of the mapping qi — ► Xi, so that the reconstruction of 
space coordinates in these local approximations is not possible (cf. subsection 2.6). 

In the basis where Jij is diagonal, 

Jij = (1 +Wi)Sij, (9) 
Raychaudhuri's equation (7) is written as 

Y , Wi x = -4irGa 4 p ( 1 - 1^ . (10) 

^(l + wt) V(! + «>i) (1 + W2) (1 + w 3 ) J V ; 

The local approximations discussed here are required to be exact for planar, spherical and cylindrical 
symmetries. In the spherical case we have w\ = wi = for a cylindrical perturbation w\ = W2, 
and W3 = 0; whereas for planar symmetry ^2 = 1^3 = 0. In these three cases, as we have only one 
independent eigenvalue u>i of the deformation tensor, this equation can be solved for toj. 



2.1. Zel'dovich Approximation 

The Zel'dovich approximation (Zel'dovich 1970), hereafter ZA, can be viewed as a solution of 
the linearized form of equation (10): 

3 3 
^2 w i = 4irGa 4 p^2wi. (11) 
i=l i=l 

Zel'dovich used the solution of the linearized equations (1 to 3) x« = qi — D (t) (q) and ex- 
trapolated it into the nonlinear regime. The eigenvalues of the deformation tensor are thus given 
by 

Wi = -D (t) A° (q) , (12) 

where A° are the eigenvalues of d^i/dx^ . Substituting this expression into the eq. (11) we find two 
solutions for D, known as the growing and decaying modes. For an EdS universe we have 

wf = - + A°a and w~ = --A°a" 3 / 2 . (13) 

Since the decaying mode becomes negligible very quickly, only the growing mode is relevant for our 
discussion. The initial conditions are specified in terms of the A?, which are functions of the initial 
positions q. The principal axes of d^i/dx^ are generally different for each point. We will denote 
the linear growing mode solution by — A^: 



A; = A°a. 



(14) 
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In the linear regime the density contrast 5 will be given by 5l = (A° + A° + A3) a = Sqo. 

The gist of the Zel'dovich approximation is that the linearized trajectories can lead to nonlinear 
density perturbations. Analogous ideas have been applied in many approximations. An example is 
the higher-order Lagrangian expansions, where the perturbed quantity is the displacement field. In 
an EdS universe the solution may be written in the form X{ = qi + Xm=i o, n ^[ n ^ (q) . The first order 

(1) . (n) 

solution is the Zel'dovich approximation. The determination of the higher order follows 
from the lower order ones through the solution of Poisson equations. The second order solution 
is known as Post-Zel'dovich approximation (Moutarde at al. 1991; Buchert 1992; Lachieze-Rey 
1993), and the third order is called Post-post-Zel'dovich (Juszkiewicz, Bouchet, & Colombi 1993; 
Buchert 1994). 

The Zel'dovich approximation is widely used for the weakly non-linear regime, and for gener- 
ating initial conditions for numerical simulations. It gives the exact solution for the case of planar 
symmetry. 



2.2. Modified Zel'dovich Approximation 

In the ZA the time factor in eq. (12) is independent of the initial conditions, and it is valid 
only for the linearized limit in Wi (eq. 11). Reisenegger & Miralda-Escude (1995) have proposed a 
generalization of the Zel'dovich approximation where D may depend on the position through the 
initial conditions A°. The Ansatz u>i = —D (r, A°) \® (q) is substituted in equation (10) to give 

d 2 D ^_ i r h D- m D 2 + m P* 

^ = ^ Gpa m - 2 m D + 3 m D* ' (15) 

where 771 = A? + A° + A3, r} 2 = A^ + + A2A3 et) 3 = A^A^Ag. This equation, which must be 
solved numerically, determines completely the function D(t, A°) and defines the Modified Zel'dovich 
Approximation, hereafter MZA. It is exact for spherical, planar and cylindrical symmetries. How- 
ever, for underdense regions (8q < 0), the MZA may not work, as pointed out by Reisenegger & 
Miralda-Escude (1995). This is due to the fact that, when not all the three eigenvalues Aj have the 
same sign, the denominator in the right hand side of eq. (15) will eventually vanish. Thus MZA 
cannot be used with this kind of initial conditions. 



2.3. Deformation Tensor Approximation 

In the two local approximations discussed above, the time dependence of the three Wi is the 
same and it can be completely determined from equation (10). The next two approximations will 
provide an equation for each of the three Wi and an analytical solution of Wi in terms of the linear 
solution Aj (eq. 14). Due to the symmetry among the axes, both the equation for Wi and the 
explicit solution in terms of Aj should be invariant under any exchange of indices, k). 
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Equation (10) may be written in the form 

3 



E 

i=i L 



(1 + Wj +Wk + WjWk) Wi - AitGa A p ( 1 + 



Wj + W k WjWk 



W; 



= 0, 



(16) 



where (i,j,k) is a permutation of (1,2,3). Audit & Alimi (1996) have, as an Ansatz, split this 
equation into three equations for each Wf. 



(1 + tfj + iy fc + T^-Wfc) ^ = 47rGa 4 p 11 + h j ^. (17) 

This equation defines the deformation tensor approximation, hereafter DTA. Another motivation 
for the above equation is that it is exact for planar, spherical and cylindrical perturbations. We 
have thus a set of local equations that allows to determine each Wi completely. Of course this 
spliting of equation (16) is not unique and we could add more local terms in equation (17) which 
would obey the symmetry requirement. 



2.4. Complete Zel'dovich Approximation 

The Complete Zel'dovich Approximation, CZA (Betancort-Rijo & Lopez-Corredoira 2000) as- 
sumes that the Wi can be expanded in terms of the linear solution Aj (eq. 14). To satisfy the 
symmetries required, the power series must have the following expression: 



r i (X i ,X j ,X k ) = l+ £ Cl^ + XkYiXj-X^Xr, (18) 

l,m,n=0 

where 

Wi = -AiT-j, (19) 

and C? „ are the coefficients of the p-th order terms, with p = I + 2n + m. The Zel'dovich 
approximation corresponds to rj = 1. The second order term 



^ (2) = -A, n (A, + A fc ) (20) 



coincides with that of the DTA. 



For planar configurations one should have rj = 1, thus Cq m = 0. The other coefficients of the 
expansion are determined from equations (1) and (3) through a recursive scheme. Betancort-Rijo 
& Lopez-Corredoira (2000) calculated explicitly the coefficients Cf m n up to the terms of fourth 
order in A in an EdS universe. 
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When the higher-order terms become important, all of them contribute roughly the same. 
Thus Betancort-Rijo & Lopez-Corredoira have chosen to truncate the series at the fourth order 
and approximate the rest by a function R(Xi, Xj, A&). This function is parametrized in such a way 
that the result is in agreement with the exact planar, spherical and cylindrical dynamics. Their 
expression for R is: 



R(Xi,Xj,Xk) — 



Xj + X), \ / Xi + Xj + Afc 



1-9 A,- J ' " 1 



2 ) \ 1.3 

^s P (Ai + A, + A fc ) - i? sp (Ai) + i? sp ( h±^i ) ) j (21) 



where R sp is the correction term R, corresponding to the spherical symmetry. By comparing the 
numerical results for overdense perturbations with the truncated series solution, they fitted R sp as 



Rsp 



(x) = 2.58xlO-V(l-^)" 1 . (22) 



The expansion (18) up to p = 4, together with expression (21), gives nearly exact results for 
spherical and cylindrical overdense perturbations, and the exact result in the planar case. Indeed, 
the CZA predicts that a spherical perturbation with Sq = 1 will collapse at a c = 1.72, whereas the 
exact solution gives a c = 1.69. 

The CZA does not apply for perturbations with negative values of A°. For example, when all 
the three A^ are negative the volume element should expand indefinitely, hence the Aj will approach 
infinity and therefore the series expansion breaks down. It can be easily seen that, as all Cf m n are 
positive, if we truncate the series in an odd power of A, r will change sign, and the fluid element 
will eventually collapse. 



2.5. Local Tidal Approximation 

The general relativistic equations for the kinematical parameters in the projection formalism 
are very akin to the Newtonian ones (see Ellis 1973). The analog of the Poisson equation is obtained 
from the equation for the Weyl tensor. Barnes & Rowlinson (1989) pointed out that by neglecting 
the magnetic part of Weyl tensor H^, the evolution equation for the electric part becomes 
local. The dynamics of kinematical parameters is then reduced to a closed set of local equations. 
This result was first applied to structure formation by Matarrese, Pantano, & Saez (1993). Since 
the magnetic part of the Weyl tensor has no Newtonian analogue, Bertschinger & Jain (1994) 
introduced the non-magnetic approximation, by simply discarding the magnetic part H^ u in the 
equation for E^ u in the application to the Newtonian cosmology. This approximation is exact for 
spherical and planar configurations, but fails for cylindrical symmetry. Also it was not able to 
reproduce the dynamics of the collapse even for a homogeneous ellipsoid. Thus we will not consider 
this approximation further in this work. 
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Bertschinger &; Hamilton (1994) pointed out that, in a "Newtonian" limit, the role of the mag- 
netic part Hij is not altogether negligible (see also Ellis & Dunsby 1997). Within this framework, 
Hui & Bertschinger (1996) have proposed the local tidal approximation (LTA), which consists in 
discarding some terms in the evolution equation for E^, to get 



dEij 1 da ^ o . . 

3 + --j-Eij = -AirGa 3 ^, (23) 



dr a dr 

where E^ is the Newtonian limit of E^ v , which gives the tidal field: 



d 2 <p _ 1 d 2 <p _ d 2 <p 4TrGa 2 p5 
dx l dx'J 3 dx k dxk dx^dxi 



E- ■ = - 5 ■ = - y §. . (24) 



Equations (23) and (6) written in terms of E^j form a closed set of local equations. It was shown 
that the LTA is exact for spherical, planar and cylindrical symmetries. In general, it is exact 
whenever the orientation and axis ratios of the gravitational and velocity equipotentials are equal 
and constant for the mass element under consideration (Hui & Bertschinger 1996). 

It is possible to show that in the LTA, once the velocity gradient is diagonalized, it will remain 
diagonal (Hui k, Bertschinger 1996) and so will the deformation tensor. We may write equation (23) 
in terms of Wi by using the equation (24) together with equation (6). We will then have a set of three 
third-order equations for Wi which completely determines their evolution, once appropriate initial 
conditions are provided. Alternatively, equation (23) can be solved in terms of the kinematical 
parameters. In this case, the Wj are calculated by (eqs. (5) and (8)) 



^ <> U • / ) «V) (25) 



where dj are the eigenvalues of the shear a. 



2.6. General Features 

The local approximations discussed above are either a system of ordinary differential equations 
or explicit expressions in terms of the linear solution. In these approximations each point evolves 
independently of the others. The influence of the other fluid elements enters only through the 
initial conditions. They give the time evolution of the deformation tensor, and thus the kinematical 
parameters for each volume element. 

These local approximations are exact under some geometrical symmetries. In particular, they 
are exact whenever w\ = W2 = w^; or w\ = W2, with ws = 0; or W2 = W3 = are satisfied locally. 
They are nonperturbative, i.e., valid, in principle, for any 5 or Aj. 
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The second order solution of the CZA and DTA (eq. 20) is in agreement with the result 
of second order Lagrangian perturbation theory for the density contrast (Sahni & Coles 1995, 
Betancort-Rijo & Lopez-Corredoira 2000). The ZA fails at second order. 

The local approximations are not appropriate for recovering the positions. To see this, let us 
consider an initial configuration such that the deformation tensor is diagonal at every point 

J ij = ^ = (l + w i )S ij . (26) 

If this holds initially, it will be valid throughout the evolution, according to the local approximations. 
In this case x« would be given by 

Xi = J (1 + Wi) d qi . (27) 

If Wi had an explicit dependence on qj or nondiagonal terms would arise in equation (26); hence 
Wi must be a function of qi only. Consequently, for this particular choice of initial conditions, each 
A? must depend only on the coordinate qi. However, as the Wi evolve, they will in general depend 
on the three A? and, ultimately, on the three coordinates: Wi = Wi (A? (<?i) , A° (52) , A3 (93) , r). 
Thus equation (26) can no longer be satisfied. This shows that the local approximations in general 
violate the integrability of the deformation tensor. In other words, we cannot recover the actual 
positions in the local approximations (except when Wi is independent of the initial position). 

Another way of seeing that the integrability is violated is as follows: If it were possible to 
reconstruct x\ from J±j, for example, then J\j ought to be a gradient field in lagrangian space. 
Therefore its curl should vanish. In the particular case of equation (26), this implies that dw\/dq 2 = 
and dwi/dq 3 = 0. These conditions are satisfied, in general, only by the ZA for which Wi = —a\f 
(provided that A° = \®(qi)). Thus, the only approximation which always permits the direct 
computation of the positions is the Zel'dovich approximation. In spite of the non-integrability, 
these methods offer an approximate solution for the deformation tensor, allowing to calculate local 
quantities, such as the kinematical parameters. 

If any eigenvalues of the deformation tensor approach —1, the density contrast S will diverge. 
Since they are functions only of a and the initial conditions, A°, we can expand them near the 
collapse time 



dwi 

Wi = -1 - — — 

da 



(a c - a) + ■ ■ ■ . (28) 



Therefore, the density contrast S behaves, at the collapse time, as 

5 oc (o c - a)" 7 , (29) 



where 7 is the dimensionality of the collapse (7 = 1 for the collapse in only one axis, 7 = 2 for the 
collapse in two axes simultaneously, and 7 = 3 for the collapse in three axes). On the other hand, 
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for the expansion 9, we have from equation (2) 



(30) 
(a c - a) 

for a — ► a c . The asymptotic behavior of the other kinematical parameters can also be determined 
in a similar fashion. 



3. APPLICATIONS 

In order to compare the performance of these approximation schemes, we apply them to some 
specific situations in the following subsections. 



3.1. The homogeneous ellipsoid 

An initially homogeneous ellipsoid in an expanding universe develops in such a way that the 
homogeneity is almost preserved during all the evolution. Therefore, the homogeneously collaps- 
ing ellipsoid model (HCE) is considered to be very accurate (Eisenstein & Loeb 1995; Hui & 
Bertschinger 1996). The result of local approximations have been compared to this model. Such a 
comparison is useful since it offers the possibility of checking these approximations in a less symmet- 
rical situation (Hui & Bertschinger 1996; Audit & Alimi 1996; Betancort-Rijo & Lopez-Corredoira 
2000). It is worthwhile to compare theses analyses including the MZA. 

The equation of motion for the HCE is given by (Icke 1973; White & Silk 1979) 

= -^al - (X 1 X 2 X 3 - Y 1 Y 2 Y 3 )C D (Y 2 ,Y J 2 ,Y 2 ), (31) 

where, as before, (i,j,k) are permutations of (1,2,3), Y; L represent the axes of the ellipsoid in 
comoving coordinates, and A« are their asymptotic values for a — > 0. The function Co is the 
degenerate case of Carlson's integral of the third kind (Carlson 1977; Press et al. 1992): 

C D {x,y,z) = - 

2 Jo (z + sy /2 y/{x + s) (y + s) 

The linear growing mode is 

Yi ^ o X t (l - i^AsXa C D (X 2 k ,Xf,X?)5 a)j . (32) 

As there is no rotation, the orientation of the principal axes does not change. Thus, the 
position of each element will be proportional to the expansion in each direction. If an element 
inside the ellipsoid has initial position qi, then its coordinates at a later time will be given by 

Y 

Xi = -z^qi- (33) 
A,; 
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With this expression we may compute the kinematical parameters which will not depend on the 
position. The same holds for the tidal field. Hence we can compute the evolution of a fluid element 
according to the local approximations, and compare with the evolution of Eij, <Jij, 5 and 9 as 
derived from the ellipsoid solution with the same initial conditions. 

From expression (33) we can see that the deformation tensor does not depend on % inside the 
ellipsoid: 



As discussed in the previous section, the positions of the fluid elements may not necessarily be 
recoverable in the local approximations. However, the choice of the same Wi for any fluid element 
(wi independent of qi) is consistent with the HCE. In this case, we may recover the positions from 
the wi as 



In Figure 1 we compare the time evolution of the axes Ri = aYi of an ellipsoid in the five 
approximations discussed in this paper. Here, the initial values Xi of the axes were arbitrarily 
chosen to be 1 : 1.25 : 1.5 with <5o = 1. The general conclusion does not depend substantially on 
the choice of these values, as will be seen in the next section. We see that the results of these 
approximations, except for the ZA, are very close to the one given by the homogeneous ellipsoid 
model. The ZA overestimates the collapse time, showing that a simple extrapolation of the linear 
trajectories underestimates the nonlinear effects. The common feature we observe in the local 
approximations is that the collapse occurs a little bit earlier in the directions of the two initially 
larger axes than the HCE case, whereas the collapse in the direction of the shortest axis is slightly 
delayed compared to the HCE. In other words, in the local approximations, the tidal forces are 
reduced compared with the HCE model. 

Concerning the collapse time a c , all these approximations give very similar results as shown in 
Table 1. The differences are less than 5%. The MZA gives the closest value to that of the HCE. 
Considering, however, that the HCE itself neglects the effect of the interaction of the background 
with the ellipsoid, this will not necessarily indicate that the MZA has the better performance among 
the other local approximations. In fact, for larger shear, the MZA deviates from the others as will 
be seen in the next section. 



Following Bertschinger & Jain (1994) we will parametrize the initial conditions in the following 

way 




(34) 



Xi = (1 + wi) qi. 



(35) 



3.2. 



Generic Initial Conditions 
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2 1 

-e Qi (a ) + -S , (36) 



where Qi (a) are the diagonal terms of the traceless quadrupole matrix 



Qij (a) = diag 



a + 27r\ fa — 2ir\ /a 
c( is j - ] . cos ( — j ( — 



(37) 



It is easy to show that Eo is related to the magnitude of the shear and tide, Qi (ao) gives ratios 
of the eigenvalues of and (note that in the linear regime E^ oc cxj.,-), and So the density 
contrast. The parameter eo varies from to oo, ao from to w, and So can go from — oo to +oo. 
However, it is sufficient to study the dynamics for Sq = +1 and So = — 1, as we shall see below. 

The initial perturbation So is defined as the ratio between S and the growth factor D in the 
linear regime: 

<5 = limA (38) 
a^O L) 

In an Einstein-de Sitter universe we have D = a. Thus choosing different values of S is equivalent 
to rescaling a. This is so for all the kinematical parameters. Therefore, the equations of motion in 
the local approximations are invariant under the following scaling 

So — ► cSq, £q —> ceq, and a — ► c _1 a. (39) 

Due to this invariance we can express the collapse time a c as (Audit, Teyssier, & Alimi 1997): 

, x \ j^lat (eo/S ,a ) , if So > 0, 

a c (So,e ,a ) = i (40) 

(Jd | a c (e / | o 1 ,«o) , if oo < 0, 

where af (eq/Sq, ao) = a c (±1, eo/ \ Sq\ , ao). Hence we just need to compute the two functions a+ 
and a~ which depend on Eq/So and ao only. 



Table 1. Collapse time for a 
homogeneous ellipsoid 



Approximation 


a c 


HCE 


1.569 


MZA 


1.578 


CZA 


1.582 


LTA 


1.612 


DTA 


1.633 
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In Figures 2 and 3, we plot the collapse time a c as a function of Eq arid «o for overdense and 
underdense perturbations, respectively. Since the MZA and CZA do not apply for some underdense 
regions, we have not displayed the results of these approximations in figure 3. We also show the 
signs of corresponding to the initial conditions in these two figures. 

The parameter space of initial conditions that can be spanned by an ellipsoid with any axes 
ratios is equivalent to having the three A° positive. The region corresponding to the homogeneous 
ellipsoid is limited to relatively small shear, and all the local approximations, except for the ZA, 
agree significantly well in this region. The ZA overestimates the collapse time for spherical con- 
figurations. We see that they still quite similar for overdense perturbations in general. The MZA 
substantially deviates from the others for high shear. In all the cases the shear accelerates the col- 
lapse, which is a well-known nonlinear effect. Thus the first regions to collapse are not necessarily 
those with higher density. We can also see that oblate initial configurations (for which cos ao > 0) 
collapse first. Thus planar collapse is favored by these approximations. 

For negative perturbations the difference among the approximations is enhanced. The LTA 
systematically gives slightly larger collapse times than the DTA. The collapse time given by the 
ZA is the shortest among the three. It is important to notice that in the local approximations 
underdense regions may also collapse, due to the effects of the shear. 

The relevance of the shear in the nonlinear phase of gravitational clustering is in agreement 
with iV-body simulations (Katz, Quinn, & Gelb 1993), yet it is sometimes ignored in structure 
formation studies. Any model based on the spherical collapse would miss this effect. 

An interesting aspect of the local approximations is that the collapse time has the same asymp- 
totic behavior 

C 

a c (5 ,e ,a ) ~ — (41) 

for high initial shear (eq S> So) in all the approximations, where C is a (slowly varying) function of 
ao only (see Appendix A). That is, the collapse time for large shear does not depend on 5o and is 
inversely proportional to the initial shear eq. 

3.3. The Cosmological Mass Function 

The mass function n (M) is defined such that n (M) dM gives the number density of collapsed 
dark matter clumps with masses between M and M + dM. These clumps are associated with 
proto-galactic haloes, and with galaxy groups and clusters. Comparing theoretical mass functions 
with observations provides important constraints on the cosmological parameters (Bahcall <fc Cen 
1993; Girardi et al. 1998; Rahman & Shandarin 2000) and the spectrum of primordial perturba- 
tions (Lucchin & Matarrese 1988; Ribeiro, Wuensche, & Letelier 2000). The approach of Press & 
Schechter (1974) to calculate the mass function, hereafter PS, has been extended to nonspherical 
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collapse and applied to some local approximations (Monaco 1995 for the ZA, and Audit et al. 1997 
for the DTA). Here, we extend such analysis to the LTA and compare them. 

Let F (M; ao) be the fraction of collapsed objects at a with mass higher than M; then the 
mass function is given by 

n (M) = -4 7 ^-. (42) 

The fraction F may be calculated as an integral over all the possible initial conditions weighted by 
their probabilities: 



*o Jo Jo J -i 



s (ao; Sq, e , ao) Pm (So, £q, «o) dSodsoda . 



(43) 



The function s is equal to one if an element with parameters Sq, £q, ao has already collapsed at ao, 
and is zero otherwise; Fo is a normalization factor. The collapse time a c of a fluid element with 
initial perturbations parametrized by So,£o,ao can be computed in the local approximations. As 
mentioned in subsection 2.6, the collapse is characterized by the divergence of the density, which 
is equivalent to the first axis collapse. Beyond this point the Lagrangian formalism breaks down. 
Some authors (Audit et al. 1997; Lee & Shandarin 1998; Sheth, Mo & Tormen 1999) have suggested 
other alternatives for the definition of collapse in the calculation of the mass function. Here we will 
prefer to keep the simplest assumption of first axis collapse, since it does not introduce any free 
parameter. 

What we need now is the probability distribution function Pm (So,£o,oto) for the initial con- 
ditions. Assuming Gaussian initial fluctuations, Doroshkevich (1970) derived the joint probability 
for the three eigenvalues of the deformation tensor \®, A[> and A3. Using this result, Pm {So, £o, a o) 
is given by the product of three independent probabilities for each parameter So, £0 and ao- 



PASo) = 



Px M 



50 



exp 



l(S_o 
2 V A 



3 V 2vrA 2 VA 



(a) 6XP ~2 (a) 



P a (a ) = sin (y) 



3 o • 2 f a 0\ 

2 sin — 

2 V 3 / 



(44) 

(45) 
(46) 



The variance A is related to the mass M and the power spectrum of the primordial density field 
Ofc through 



A 2 (R) = / 
Jo 



dk 
2^2 



Wl(R)k 2 a\, 



(47) 
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where M = (4ir/3)fwR 3 P and Wk (R) is the Fourier transform of a filter with width R in physical 
space. The factor fyy depends on the shape of the filter function, for a top-hat filter we have 
fw = 1, whereas for a sharp— filter fw = 9ir/2. The mass function can now be written in the 
form 

where <3? (A) = dF(A)/dA. The function <E> (A) contains all the influence of the dynamics and 
depends neither on the particular form of the power spectrum nor of the filter W; it is referred to 
as the universal mass function (Audit et al. 1997). 

We calculate the universal mass function for the ZA, LTA and DTA but not for the MZA 
and CZA since they do not apply for negative density perturbations. In appendix B we show the 
detailed calculation. 

In Figure 4, we show the mass functions for those approximations. For comparison, we also 
display in this figure the fit to A-body simulations obtained by Jenkins et al. (2000), together with 
the standard PS mass function. 

We see that the results of the DTA and LTA are very similar. Furthermore, in the high-mass 
tail (A < 0.5), they reproduce well the results of the A-body simulations. However, we can see 
that these approximations overestimate the concentration of masses near A = 1. The right-end 
tail of the distribution decays more rapidly compared to the A-body simulations. This tendency 
is still enhanced in the ZA. However, in these approximations, the position of the maximum of the 
distribution is close to that of the A-body simulations, giving a better estimate than that of the 
PS; in particular, the LTA and DTA give nearly the same value as the A-body results. 

As for the normalization factor F , there exists an extensive discussion on its origin (see, for 
example, Peacock & Heavens 1990; Bond et al. 1991; Jedamzik 1995; Yano, Nagashima, & Gouda 
1996). The normalization factors for the local approximations are close to one (1/0.92 for the DTA, 
and 1/0.89 for the LTA) whereas in the original PS derivation the normalization factor needed is 
Fq = 1/0.5. This is due to the fact that, in the spherical collapse model, only overdense regions 
collapse. 

The fact that around A = 1 the number of objects is overestimated in the local approximations 
implies, due to the normalization of the mass function, that they should provide a lower estimate 
than the A-body simulations for large enough A. There, the contribution from the low-mass objects 
is dominant; in any realistic process, they may also arise from the fragmentation of larger clusters. 
The criterion for the formation of a clump from the direct collapse of an initially perturbed region 
does not account for these complex processes of fragmentation. Therefore, the discrepancy of the 
mass function for high A might be attributed to the use of expression (43) rather than to the 
definition (42). 



4. DISCUSSION 



We have investigated local Lagrangian approximations to the nonlinear dynamics of pressure- 
less dark matter. We have selected the modified Zel'dovich approximation (MZA), the deformation 
tensor approximation (DTA), the complete Zel'dovich approximation (CZA), and the local tidal 
approximation (LTA), in addition to the original Zel'dovich approximation (ZA). These four ap- 
proximations were designed to improve the ZA, and are in fact exact for planar, spherical and 
cylindrical symmetries, whereas the ZA is only exact for the planar case. They are semi-analytic 
and easy to be implemented in any application where local quantities are involved, such as the 
calculation of the mass function in the PS approach. 

All the local approximations discussed here, except for the ZA, provide quite a similar evolution 
for an ellipsoid, reproducing the results of the homogeneous ellipsoid model. Thus, for this kind 
of positive density perturbations, these methods work fairly well. However, the MZA turns out 
to deviate substantially for large values of the shear as was shown in section (3.2), reflecting the 
fact that it does not give the correct second order solution. Furthermore, the MZA cannot deal 
with initially underdense regions that will eventually collapse. Therefore, its applicability is rather 
limited when compared to the other approximations. 

We note that the second order expansions of the CZA, LTA and DTA coincide with the second 
order Lagrangian perturbation theory, whereas those of the MZA and ZA do not. It is interesting 
to recall that the LTA and DTA have very different origins from the CZA, but still they give the 
correct second order result. 

The CZA, LTA and DTA give quite analogous results for generic initial conditions. However, 
at least in its original form, the CZA cannot be used for negative values of Aj. One possible solution 
to this problem might be achieved through an expansion such as 

_ 1 + Y2l!m,n=0 E l,m,n(^j + (Aj - Afc) 2n Af t 

n ' 1 + £f m ,n=0 SUA- + " \ k ? n \? ' 

The coefficients E and D should be appropriately chosen to adjust the asymptotic behavior; in par- 
ticular, we could use the numerical solution for underdense cylindrical and spherical perturbations 
to fit some of these coefficients, as done for the overdense case in the CZA. Besides, to agree with 
the perturbative solution (18) we should have 

For higher orders, the determination of these coefficients is rather complicated. Further investiga- 
tions on this possibility should be pursued. 

Concerning the mass function, it is found that the LTA and DTA give an accurate result for 
large masses as compared to the A^-body simulations. The position of the peak is also in good 
agreement, whereas its amplitude is overestimated by a factor 2 compared to the iV-body results. 
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Since the mass function is normalized to unity, this means that the local approximations, together 
with the PS formalism, underestimates the density of low-mass clusters. However, this might be a 
consequence of the criterion for the formation of a collapsed object based only on the collapse time. 

It is interesting to notice that the collapse time, as a function of «o an d £q, has an approximate 
scaling property (see eq. Al), which is very precise for £q/8q 3> 1. We conclude that this may be a 
general feature of gravitational collapse in local approximations, whose validity is worth checking 
in a more general setting. 

While there is still not a clear theoretical understanding or support to the local approxima- 
tions, they proved to be very accurate in the situations investigated here. They reproduce some 
well known features of nonspherical collapse, such as the possibility of collapse of some initially 
underdense regions, and the fact that the shear accelerates the collapse (see Sahni & Coles 1995). 
The main limitation of these approximations is that they only provide information about the in- 
ternal state of a given mass element, but do not determine its position. Even so, their simplicity is 
highly expedient for practical applications, such as the calculation of nonlinear corrections to the 
microwave background anisotropics, and the Gunn-Peterson effect. In particular, they are suitable 
for obtaining statistical properties of the present fields as a function of the primordial ones, as in 
the case of the mass function. Further studies on the validity of these approximations, based on a 
comparison to iV-body simulations, are required. Such a comparison would allow to fully test the 
approximations described in this paper and more generically, the locality hypothesis. If they still 
provide accurate results in this case, the local approximations could represent good alternatives to 
the computer simulations, taking much less computational time, allowing thereby a larger scan- 
ning of initial conditions. They could give complementary information to the iV-body simulations 
and would provide a better physical understanding of the nonlinear dynamics of self gravitating 
systems. 

Most of the results of this paper may be extended to more general backgrounds. The influence 
of any smooth component only alters the behavior of a(t), and of the growing mode growth factor 
D (t), that will not be equal any more. It would be interesting to study the relativistic analogue of 
the LTA. Another interesting extension would be to include vorticity in the local approximations 
as was done for the ZA in (Buchert 1992; Barrow & Saich 1993). We could use them to test the 
effects of a possible primeval vorticity on large scales (Li Xin-Li 1998). 
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A. FITTING FORMULAE FOR THE COLLAPSE TIME 

In order to avoid repeated numerical integrations of differential equations in the LTA and DTA, 
we have parametrized the collapse time as a function of initial conditions. For both of these local 
approximations, the following scaling property is approximately satisfied: 

af (x, ao) ^ # ± (/) , with / = xg ± (a ) , (Al) 

where H ± and g ± are functions to be fitted for each approximation, and x = Eq/ \8q\. This relation 
becomes more accurate for increasing x. In any case the error of the fit is less than a few percent. 

For g (a) we found that a kind of truncated Fourier series can be used to a very good approx- 
imation: 

g (a) = c\ cos ^— ^ + c 2 cos a + C3 cos + C4 cos (^ a ) + 

c 5 cos ( J + (1 - ci - c 2 - c 3 - c 4 - c 5 ) cos (3a) . (A2) 



x 2 , 

The values of the parameters for each case are shown in Table 2. 

For the dependence of a c on x for fixed ao, we fitted the function / = H^ 1 rather than H 
itself, because this is the quantity we will need to compute the mass function. Taking into account 
the boundary value and asymptotic behavior, we parametrize / by 

, diz (l + d 2 z + d A z 2 ) 

f = , , — 1 ox , with z = d - a c , A3 

(1 + d 3 z + d 5 z 2 ) a c 



for overdense regions (Sq > 0), and 



/" = * + ^ _ 2 ^, (A4) 

1 + d 3 a c + d 5 a c 



Table 2. Parameters of p(a) fitted for the 
LTA and DTA 



9 


Cl 


C2 


C3 


c 4 


C5 


9lta 


1.546 


-1.015 


0.786 


-0.462 


0.182 


9lta 


1.461 


-0.767 


0.473 


-0.231 


0.079 


9dta 


1.505 


-0.842 


0.508 


-0.228 


0.066 


9dta 


1.497 


-0.836 


0.513 


-0.234 


0.068 
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for underdense regions (do < 0). The parameters for the LTA and DTA are given in Table 3. 
Notice that for high shear (x S> 1) we have a c Cl such that: 

+ did (l + d 2 d + cMo) 1 ... 
° c ™ (1 + d 3 d + d 5 d$) xg+(a )' 

and 

d\d 4 1 . . 

a c — ► — — — -, A6 

x^co d 5 xg (a ) 

In the Zel'dovich approximation we have an analytical expression for a c . As the collapse occurs 
when the greatest Aj reaches the value 1, using the parametrization (36) we get: 

= 5 + 2e cos(a /3)- (A7) 
In this case we can clearly see the features of a c : 

ZA 1 S 

a c (6 ,e ,a ) ~- 



' 1 "" l |<5 | ±1 + 2 ( £o / |«5 |) cos (cio/3) 
1 



j^a^ A (±l,e /\5o\,a ). (A8) 
|oo I 

Note that the property (Al) is satisfied exactly for the ZA. 

B. CALCULATION OF THE MASS FUNCTION 

With the integral (43) we may write the universal mass function <I> (A) = dF/dA in the form 
1 d f n f°° f^ 00 

$ (A) = Fh dAJo Jo J s (v^X^ao)P(v,X,ao)dvd X dao, (Bl) 



Table 3. Parameters of / fitted for the LTA and DTA 



/ 


do 


di 


d 2 


d 3 


C?4 


d 5 


flTA 


1.686 


8.469 


19.88 


78.49 


17.03 


163.2 


f DTA 


1.686 


14.14 


13.34 


87.68 


8.676 


163.8 


f LTA 


0.591 


1.064 


0.678 


-1.335 


16.054 


8.613 


f DTA 


0.495 


0.942 


0.322 


-2.083 


25.718 


12.961 
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where v = So /A and \ = £o/A. With these new variables the dependence on A will be present 
only in the function s, which may be written as 



s = (1 - a c (5q,£q, a )) 



(B2) 



Note that in the case of the Press & Schechter original approach this function is given by: s 
(5q — S c ) , where S c = 1.686 is the value at which a spherical perturbation collapses at a = 1. 



To calculate ds/dA one uses the relation (40), obtaining: 



ds 
dA 



5 D (1 - a c (iyA,xA,a )) 



MA 2 



(B3) 



where So is the Dirac delta function. Therefore we may eliminate one of the integrals in expression 
(Bl), with the mass function being calculated over the surface a c (vA, %A, ao) = I^Ap 1 a^r (x/V, olq) 
= 1. For this sake we need to write one of the three variables in terms of the others on this surface. 

Let us assume the we have x as a functions of vA and a: X = Xa (^A, a), where the subscript 
a indicates that x is calculated over the surface a c = 1. The integral in x m equation (Bl) is thus 
eliminated using the relation 



S D (1 - a c (vA, xA, ao)) = Sd(x~ Xa) 
The mass function (Bl) is now given by: 



dx 



$(A) 



^ rir r+oo ^ 

To Jo i-oo A 



da f 



dx 



Pv (v) P x (xa) Pa (ao) dudao, 



(B4) 



(B5) 



where da c /dx is evaluated in Xa ("A, a) . We can simplify this expression further if the collapse time 
a c (x, a) is only a function of the product xg (a) , which we have seen is an excellent approximation 
for the local approximations studied here (see appendix A). Using the property (Al) we have 



Xa 



g(uo) 



H O a I) > 



and 



da c 




1 dHg(a ) 


dx 




vA df v 



(B6) 



(B7) 



where the superscript + is implied for positive v, and the — for negative v. Replacing these results 
in expression (B5) we get finally: 



d>(A) 



1 



f dH 



1 



g(ao) 



p v (v) P x (Xa) P a (ao) da du, 



(B8) 



where Xa = W\ fa/g (ao), with f a (uA) = H 1 (|^A|) and dH/df, is calculated in f a . This is why 
we have chosen to fit the function H~ l , instead of its inverse. 
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As P (v) and dH/df are independent of oq, we integrate first in this variable: 



r i 

h(u,A):= / - T -— 
Jo ff 5 ("o) 

The universal mass function will now be given by 

i r+oo 

*{A) = —N / h(u,A)u 6 f 4 a 

^0 J-oo 

where N = 50^5/6 it is the product of the normalizations for P„ and P x (equations (44) and (45)). 
Note that the above integral is limited for positive values of v, as H^ 1 (|^A|) = for u/S. > /q 
(fo = 1.686 for the DTA and LTA, and f = 3, for the ZA). Although our fitting formulas (sec. 
A) become less accurate for v — > — oo and v — ► /o/A, the mass function is not affected, since the 
integrand in (BIO) goes to zero in these regions. Here it is clear that the underdense regions do 
contribute to the mass function, as pointed out by Audit et al. (1997). 

As an example, let us consider the Zel'dovich approximation. In this case we have 

H± (/) = ±TT2/' and 9 {a) = cos (f ) • (B11) 

Using the following transformation of variables 



g(ao) 



P a (a ) da . 



(B9) 




exp 



v 
~2 



dv. 



(BIO) 



sin (^j — > P a da = x 



3 n 

2r 2 

2 ZX 



dx 



we find an analytical expression for the integral (B9): 



1 



■ exp 



\<yf? 



1 



25 (ufY 



exp 



15 f n 2 



_9_ 
10 



where f± = (3/ A) T 1) /2. The mass function will be given by: 



(B12) 



(B13) 



. . 1 15^5 f +QO 1 / v 2 \ 

*(A) = ^— y^ ^ exp v"Tj exp 



F 8vr 
12 



5 (3 - z/A) 2 
'8 A 2 



-A z exp 



15(3-z>A) s 
¥ A 2 



l)+^(3-^A) 2 



(B14) 
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Fig. 1. — Evolution of the three axes Ri of an ellipsoid according to the homogeneously collapsing 
ellipsoid model (HCE, solid curve), and the five approximations considered in the text: Zel'dovich 
(ZA, dotted curve), Modified Zel'dovich (MZA, short dot), Deformation Tensor (DTA, dash-dot), 
Complete Zel'dovich (CZA, short dash) and Local Tidal (LTA, long dash). The initial axes ratios 
are 1 : 1.25 : 1.5 and the density contrast linearly extrapolated to a = 1 is 5q = 1. The ZA 
overestimates the collapse time, whereas all other approximations are close to the HCE. 
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Fig. 2. — (a) The collapse time as a function of the initial conditions for overdense perturbations 
with 5q = +1. The contours of constant collapse time, expressed by the scale factor a c , are displayed 
for the ZA. The light (heavy) contours are spaced by 0.1 (0.5) in a c , with the outermost contour 
being a c = 0.4 and the central velue a c = 3. (b) The same as (a) except that the MZA is used. 
The innermost contour is a c = 1.6. 
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Fig. 2. — (c) The same as (a) except that the DTA is used. The innermost contour is a c = 1.6. (d) 
The same as (a) except that the CZA is used. The innermost contour is a c = 1.7. (e) The same 
as (a) except that the LTA is used. The innermost contour is a c = 1.6. (f) The signs of A^ in each 
region of the parameter space of initial conditions, with 5q = +1. The inner region corresponds to 
the values of \® that can be spanned by a homogeneous ellipsoid with any axes ratios. Spherical, 
planar, and cylindrically symmetric perturbations are marked with a cross, a plus sign, and a star 
respectively. The upper part of these graphs, which correspond to < ao < vr, have A3 > A° > A?. 
The lower one covers the same values of A°, through a permutation of the indices. 
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Fig. 3. — (a) The collapse time as a function of the initial conditions for underdense perturbations 
with So = —1. The contours of constant collapse time, expressed by a~ , are displayed for the 
ZA. The light (heavy) contours are spaced by 0.1 (0.5) in aj , with the innermost contour being 
a^T 1 = 0. Initial perturbations in the central region do not collapse, (b) The same as (a) except 
that the DTA is used, (c) The same as (a) except that the LTA is used, (d) The signs of A- 1 
corresponding to each region of the parameter space of initial conditions, with 5o = —1. Spherical, 
planar, and cylindrically symmetric underdense perturbations are marked with a cross, a plus sign, 
and a star respectively. The upper part of these graphs, which correspond to < ao < vr, have 
A3 > A<] > A°- The lower one covers the same values of A- 1 , through a permutation of the indices. 
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Fig. 4. — (a) The universal mass function calculated for the ZA (dotted curve), the DTA (dash- 
dotted curve) and the LTA (dashed curve). For comparison, we display in this figure the fit to 
iV-body simulations (line+cross) obtained by Jenkins et al. (2000), together with the standard PS 
mass function (solid curve), (b) The same as (a) but now a logarithmic scale is used in the y axis. 
Here the x axis scale is limited to the ranee 0.332 < A < 3.32 covered bv the iV-bodv simulations. 



